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ABSTRACT 

A pseudo-time method is introduced to integrate the compressible Navier- 
Stokes equations to a steady state. This method is a generalization of a 
method used by Crocco and also by Allen and Cheng. We show that for a simple 
heat equation that this is just a renormalization of the time. For a 
convection-diffusion equation the renormalization is dependent only on the 
viscous terms. We implement the method for the Navier-Stokes equations using 
a Runge-Kutta type algorithm. This enables the time step to be chosen based 
on the inviscid model only. We also discuss the use of residual smoothing 
when viscous terms are present. 
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1 . INTRODUCTION 


The solution of the compressible Navier-Stokes equations for flow about 

two- and three-dimensional complex aerodynamic configurations is still a time 

consuming problem on today's supercomputers. The resolution of the boundary 

layers requires the use of very fine meshes in the neighborhood of solid 

bodies. For a typical viscous flow the mesh can be several orders of 

magnitude finer (depending on the Reynolds number) than that required for an 

invlscid calculation. As an example, using a C-type mesh about an NACA 0012 

airfoil, a typical mesh spacing near the body in the normal direction for an 

-2 

inviscid calculation is 1 * 10 chords. For a laminar viscous calculation 

3 -4 

with Re = 5 x 10 , this minimum cell height would be about 6 x 10 

chords. For a turbulent calculation using an algebraic turbulence model and 

with Re = 3 x 10^, the minimum cell height would be about 8 x 10 ^ 

chords. In all cases a typical chordwise spacing at the midsection of the 

-2 

airfoil is about 5 x 10 chords. 

Using an explicit method this fine mesh reduces the time step, due to 
stability requirements, that can be used. The time step restriction is caused 
by two factors. One contribution is due to the effect of the finer mesh on 
the inviscid portion of the calculation. When using an explicit method this 
reduction of the time step cannot be avoided without using a coarser mesh. It 
follows strictly from the need to include the entire domain of dependency in 
the numerical algorithm. Use of a local time step allows faster convergence 
to a steady state, but it does not remove the requirement to satisfy the 
convection stability condition in a local sense. A second difficulty is 
caused by the viscous terms. For an explicit method the time step is now 
dependent on the square of the mesh size rather than just the mesh size as 



occurs for inviscid flow. Thus, even for a high Reynolds number flow the 
viscous time step will dominate when the mesh is sufficiently fine. In all 
these cases the use of an implicit scheme will alleviate the difficulties. In 
some ADI methods the Jacobian of the viscous terms is not used in the implicit 
portion of the code in order to improve the speed of the calculation [7], We 
thus conclude that for both explicit and many implicit codes it is 
advantageous to account for the dependence of the time step on the viscous 
terms. 

In this study we shall only discuss steady state problems which are solved 
by a pseudo time-dependent method. Hence, we can change all time derivatives 
as long as the steady state solution is not affected. One common device is to 
use a different time step in each zone. It is easier to calculate this local 
time step based on the inviscid equations. This provides an additional reason 
to eliminate the dependence of the time step on the viscous terms. 

In this study we shall analyze a method used by Crocco [4] and also by 
Allen and Cheng [2]. They claim that the new scheme is unconditionally stable 
for a simple diffusion equation. We will show that in effect the scheme is a 
standard Euler forward-in-time central-in-space scheme. The time is 
artificially slowed down so as to satisfy the stability criterion. We then 
extend this scheme to the compressible Navier-Stokes equations using a Runge- 
Kutta scheme [9]. This modification enables us to choose our time step based 
on the inviscid equations. The modification automatically reduces the local 
time step in regions where the viscous time step is of importance. This 
enables us to use the inviscid time step in the far field while automatically 
accounting for viscous effects in the boundary layer. We will also look at 
residual smoothing for the heat equation. 
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II. SCALAR EQUATION 

In this section we analyze and extend a scheme for the Navier-Stokes 
equations proposed by Crocco [4] and Allen-Cheng [2]. This scheme was also 
analyzed by Peyret and Viviand [6] and Roache [8], and we will extend their 
analysis. 

We first consider the heat equation 
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Thus, for this model problem the Crocco scheme is identical with the Euler 
scheme (2) with an artificial time step At g given by 
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Thus, the unconditional stability is achieved by slowing down the time 

2 

process. Note that as At ■*■ », At (Ax) /2e, i.e., the stability limit for 

e 

the Euler method. So choosing a large time step for (4) is equivalent to 
choosing At g at the stability limit for (2), and we have merely scaled the 
time. This can also be derived from the modified equation given in [6] . If 
e or Ax is not constant, this also introduces a local time step. 

We next consider the convection-diffusion equation 


w = aw + ew . 

t X XX 


( 8 ) 


The Crocco scheme now becomes 
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with At given by (6). Thus, again this is equivalent to the Euler scheme 
with a time scaling that depends only on the viscous terms. Allen and Cheng 
utilized this scheme within a time-marching scheme proposed by Brailovskaya 
[3]. We generalize this by considering a general N-stage Runge-Kutta scheme. 

Consider the two-dimensional equation 


w^=Hw+ew + e„ w (11) 

t 1 xx 2 yy 


where Hw describes the hyperbolic or first-order terms. In [9] we describe 
a Runge-Kutta scheme where the viscous terms are frozen for all the stages. 
This is similar in philosophy to the Brailovskaya scheme. Using the Crocco 
formulation the (K + l)-st stage becomes 
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This reduces to a Runge-Kutta scheme 
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where Hq, are the approximations to the hyperbolic and parabolic parts 

respectively and 
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We slightly generalize (14) by redefining At g by 


_1_ 

At 
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( 15 ) 


where < is a constant that we can choose. The form of (15) no longer 

follows directly from the Crocco formulation. Instead k will be chosen 

based on a stability analysis. 

We choose At in (12) or (15) based on the hyperbolic (inviscid) 

stability condition. We then find At g from (15) and advance to stage 
(K + 1) using the Runge-Kutta scheme (13). 

The constant k in (15) can be chosen so that we recover the parabolic 
stability limitation when H D = 0. The exact value of k depends on the 
coefficients u^ in the Runge-Kutta formula. In order to see this more 

clearly we revert to the one-dimensional convection-diffusion equation (8). 
We replace all space derivatives by second-order central differences while the 
time derivative is kept continuous. We therefore have 


r n n 'i 
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(16) 


We Fourier transform (16) to get 


w^ *= Xw 


with 


A(5) = - ■ — (1 - cos g) + sin ? ° < ( < 2tt . 

(Ax) 2 AX 


(17) 


(18) 


A Runge-Kutta scheme for (16) or (17) is stable whenever z(g) = X(5)At g lies 
within the stability domain that depends on a i>***> a N for all 0 £ < 2tt. 





We consider the stability domain for the four-step scheme with = 1/4, 

= 1/3, = 1/2, = 1. This scheme has a stability condition along the 

imaginary axis of max|z| 2/2, i.e., for a hyperbolic problem (e = 0) 
aAt _ ^ 

- — — - 2/2. Along the negative real axis the stability condition is 

2eAt 

|z| < 2.8 and for a parabolic problem (a = 0) < 2.8. Hence for this 

(Ax) Z ~ 

case we would choose k in (15) as k = 1.4. We define the cell Reynolds 
number as 



(19) 


The previous analysis shows that the Runge-Kutta scheme is stable for R h = 0 
and R^ = We do not have any proof that the scheme is stable for all R h . 


III. NAVIER-STOKES EQUATIONS 

We now discuss the implementation of these ideas to the two-dimensional, 
compressible, Navier-Stokes equations. The extension to three dimensions is 
straightforward. We first consider the conservation form in Cartesian 
coordinates. We express the equations in the following form 
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and Hj denote first derivative terms (including the artificial viscosity and 
also the viscous dissipation function). The coefficients of viscosity (p 
and X), Y the specific heat ratio, and the Prandtl number Pr are all 
assumed (for the analysis) to be locally constant. 

In deriving our results we shall ignore all cross derivatives (see, e.g., 
[1], [2]). Based on our previous analysis we add the following terms to the 
standard Runge-Kutta scheme. 
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where Aw = w n+ ^ - w 11 and denote the usual space derivative terms. 

For simplicity we have chosen k = 1, and a denotes the constant in the 
Runge-Kutta scheme (28). Thus the density equation is unchanged. The second 
and third equations can be solved directly for A(pu), A(pv). Once A(pu), 
A(pv) are known the last equation can be solved for A(pE). As before these 
corrections imply an effective time step which automatically accounts for the 
viscous time step. In this case the effective time step differs for each 
equation. 

We finally consider the Navier-Stokes equation in body fitted 

coordinates. This can be done either in a finite volume scheme or by using 
transformations. The result is the same in either case [9], and so we shall 
use a transformation for ease of presentation. Let C = 5(x,y), n = n(x,y) 
be the body fitted coordinates. We choose the coordinate scaling so that 
AC = An = 1. The Navier-Stokes equations (20) now become 
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where h\ are first derivative terras and we have ignored all second cross 
derivative terms. As before this generates an appropriate correction term to 
the Runge-Kutta scheme. Equation (21) is now replaced by 
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where represents the standard finite difference terms. 

As before the density equation is unchanged by the viscous correction. 
Now, however, the two momentum equations are coupled together, unless the 
coordinate system is orthogonal. As we have two equations for A(pu) and 
A(pv), and we can easily solve these. To simplify the notation we define 


z = 1 + 


2aAt 
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A p = K, 


A(pu) = 


K 2 z 4 - h Z 2 


A(pv) 


K 3 Z 1 “ K 2 Z 2 


As before given A(pu) and A(pv) we can solve for A(pE) directly from the 
energy equation in (23). We also note that if one uses the thin layer 
approximation (dropping all second £ derivatives and cross derivatives in 
(22)) then these terms simplify slightly. In this case Ap, Apu, Apv are 
still given by (25) with 

z l - 1 + — — [U + 2p)n x + un y J 
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IV. RESIDUAL SMOOTHING 

As an alternative method of reducing the effect of the parabolic terms on 
the stability of the scheme we consider residual smoothing. With this 
technique one post-processes an explicit method with an implicit method. In 
practice one post-processes each equation separately and each direction 
separately so that only scalar tridiagonal matrices need be inverted. When 
using a multistage Runge-Kutta method, one can apply the residual smoothing 
after each stage, or at the end of the entire process, or any intermediate 
permutation. 

In [10] it is shown that one can construct such a scheme for a hyperbolic 
equation so that the total method is unconditionally stable. It is further 
shown in [10] that it is not efficient to use a very large At even ignoring 
splitting errors. An optimal At is about two to three times larger than the 
explicit time step. We now consider the process for a parabolic problem in 
order to see the effect of viscous terms. 

We, therefore, consider the heat equation 


u = bu 
t xx 


We solve this equation by a k-stage Runge-Kutta scheme 


(1) n , 

u *• u + c*j AtQu 
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where are given coefficients with = 1. Q is a difference 

approximation to u^. The amplification factor corresponding to (28) is 


G - 1 + AtQ + 8 2 (At) 2 Q 2 + ••• + B k (At) k Q k (29) 


where B^ = 1 and B^ = B^_^ °i c -£+i » * “ 2,»«»,k. Q is the Fourier 

transform of Q. Hence, for second-order central differencing 


Q = " 


4b sin (9/2) 
(Ax) 2 


(30) 


Residual smoothing consists of updating a stage (£) by 


/ * _ v. (£) (£) n 

(1 - aD 2 )Au = u - u 


(31) 


where D 2 is again a second-order central difference approximation to u xx , 
i.e., D 2 + (1,-2,1). We now consider two possibilities. In the first we 
apply (31) only after the final stage. Then the new amplification factor is 


GjO) = 1 + 


Bj AtQ + B 2 (At) 2 Q 2 + ••• + B k (At) k Q k 
1 + 4a sin 2 (6/2) 


(32) 


The second case we consider is applying (31) after every stage. The resultant 
amplification factor is 


G 2 (9) = 1 + B 1 AtR + B 2 (At) 2 R 2 + ••• + S k (At) k R k 


(33) 


with 


R = 


1 + 4a sin (9/2) 


(34) 
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We now investigate the possibility that either of these schemes is 
unconditionally stable. To investigate this we need only consider At 
sufficiently large. We thus consider At * « with a -*■ ®. Then (32) becomes 

(-l) k B k [^ ? si„ 2 (i)] k 

G. (e) — ► 1 + ( A x) . (35) 

1 + 4a sin (0/2) 

We thus see that for k even, G^(0) > 1 and so (28) - (31) cannot be stable 
for At large. For k «* 1 the scheme is identical with backward Euler for a 
scalar one-dimensional equation and, hence, unconditionally stable. For the 
second case we see that (33) has the same form as a standard Runge-Kutta 

A A 

method with Q replaced by R, (34). Hence, it follows that the scheme is 

A 

stable whenever AtR is within the stability region of the scheme. As 
At ®, so does a and so there is a cancellation between the numerator and 

A 

denominator; thus, AtR remains bounded as At increases. We thus conclude 
that applying the residual smoothing after each stage can make the scheme 
unconditionally stable even for a Runge-Kutta method with an even number of 
stages. 

We also see from the above argument that as At increases so must a. In 
[9], [10] we show that for a hyperbolic equation 

Uj. + au x = 0 
2 

that 0 is proportional to (aAt/Ax) . For the parabolic problem (27) it 
follows from (35) that a should be proportional to bAt/(Ax) . For the 
combined convection-diffusion equation a will be related to the sum of two 


such contributions 
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It follows from (33), (34) that if we apply residual smoothing after every 
stage then the stability polynomial has the same form as the original 

A A 

polynomial (29). The only difference is that Q is now replaced by R. From 

A A A 

(34) it follows that the ratio of Q to R is real. Hence, if Q is any 

A 

complex number then R lies along the same ray in the complex plane but with 
a different amplitude. We therefore have shown that if the original scheme 
was unstable for a given direction then residual smoothing cannot stabilize 
the scheme. Furthermore, if the original scheme was conditionally stable then 
by choosing a = o(At) sufficiently large we can make the scheme 
unconditionally stable. We have thus shown 

A 

Theorem: Let Q be the amplification factor for any approximation to the 

convection-diffusion equation and let (29) be the stability polynomial for a 


k stage 

Runge-Kutta 

scheme. 

We i 

now apply residual smoothing, (31), after 

every stage of the 

scheme. 

If 

the original scheme was unconditionally 

unstable 

then the new scheme 

is 

still unconditionally unstable. If the 

original 

scheme was 

conditionally 

stable then the scheme with residual 


smoothing can be made unconditionally stable by choosing o(At) sufficiently 
large . 

Hence, if the smoothing is applied at the end when solving a parabolic 
equation, then the scheme can be unconditionally stable only when using a 
multistage scheme with an odd number of stages. When the smoothing is done 
after each stage, the scheme can be stabilized for a large. For a system 
with a hyperbolic portion and a small parabolic contribution, e.g., high 
Reynolds number Navier-Stokes , the residual smoothing is most effective with a 
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time step about twice that of the explicit convective portion. Hence, the 
question of unconditional stability is somewhat academic. In practice [8] the 
Runge-Kutta scheme for the Navier-Stokes equations is used with four stages 
and with the residual smoothing applied after each stage. 


V. RESULTS 

In this section we present some results for viscous flow obtained using 
the analysis of Sections II and III. We used a Runge-Kutta code to solve the 
Navier-Stokes equations for two flows about an airfoil section. The details 
of this code are discussed ip [5], [9], [10], [11], In these cases we 
considered only the thin-layer form of the Navier-Stokes equations. 

For the first case we computed laminar flow over an NACA 0012 airfoil with 
a free-stream Mach number of 0.5 and a Reynolds number Re^ of 

3 

5 x 10 . The angle of attack (a) of the airfoil was zero degrees. Half- 

plane calculations were performed using a C-type grid consisting of 64 cells 

in the streamwise direction and 64 cells in the normal-like direction. The 

-4 

grid spacing at the airfoil surface was about 6 x 10 chords. The mesh 
spacing in the streamwise direction over the central part of the airfoil was 
AX = 0.05 chords. Results for this case are shown in Figures la - lc. As 
indicated in Figure lb, the flow separates at X = 0.817 chords. The size of 
the recirculation zone is displayed in Figure lc. The results are all 
independent of the time step procedure used to reach the steady state. 

In Figure Id convergence histories for this case for two calculations are 
shown. The residual displayed in this graph is the root mean square of the 
residual of the continuity equation. The calculations were started 
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impulsively by inserting the airfoil into a uniform flow and immediately 
enforcing the appropriate boundary conditions. Local time stepping and 
enthalpy damping, (see [9]) were employed in each computation; no residual 
smoothing was used. For history A the Runge-Kutta scheme with the time step 
(At) limitation determined by convection was used; this required choosing a 
CFL = 1.0. For curve B a larger Courant-Friedrichs-Lewy (CFL) number was 
used by accounting for the diffusion limit on At with the pseudo-time 
algorithm. This allowed choosing CFL = 2.5 based on an inviscid 

criterion. There is additional work with the pseudo-time scheme. 
Nevertheless, the computational time required to reach a satisfactory level of 
convergence was reduced by a factor of 1.7. 

In the second case we solved for turbulent flow over an NACA 0012 airfoil 
with M m = 0.5, Re^ = 2.89 x 10^, and a ■ 0 degrees. A 60 x 50 half-plane 
grid was used in the computations. The grid spacing at the surface was about 
8.5 x 10 5 chords. The chordwise spacing at the raidsection of the airfoil 
was about AX = 0.036 chords. Numerical results for this case are presented 
in Figures 2a and 2b. 

Figure 2c shows two convergence histories for this turbulent flow case. 
As in the laminar flow problem, the histories were obtained by computing 
without and with the effects on At due to diffusion. The pseudo-time 
algorithm was about 1.4 times faster in reaching steady state. This is close 
to the factor expected, since we were able to increase the CFL from 1.5 to 
2.7, a factor of 1.8. We do not achieve this speedup of 1.8 since there is 
some reduction of the effective time step due to the diffusion terms. 
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VI. CONCLUSIONS 

The use of the Crocco scheme for a scalar convection-diffusion equation 
introduces a scaling of the time step. This reduces the effective time step 
so that the viscous stability limit is automatically satisfied. As such the 
scheme cannot introduce any fundamental acceleration in reaching the steady 
state. The advantage of the scheme is that we do not need to explicitly 
account for the viscous time step restriction; it is done automatically. This 
can be done efficiently using Runge-Kutta type schemes. In addition, for 
variable coefficients or nonuniform meshes this introduces an effective local 
time step. 

Using this scheme for a system of equations, e.g., Navier-Stokes , has the 
additional benefit that a different scaling is chosen for each equation. Thus 
each equation has its own appropriate (viscous) time step. This is equivalent 
to using a diagonal preconditioning [10] to accelerate the equations to a 
steady state. Computations demonstrate that we can gain a factor of between 
1.5 and 2 with little programming effort. 

We further show that if one uses residual smoothing to increase the time 
step then one must also account for the viscous terms. When the smoothing is 
applied after the completion of a Runge-Kutta cycle then unconditional 
stability is possible only if an odd number of stages is used. Applying the 
smoothing after each stage allows for unconditional stability for all 
multistage schemes provided a is chosen sufficiently large. 



- 20 - 


REFERENCES 

[1] S. Abarbanel and D. Gottlieb, "Optimal time splitting for two- and 
three-dimensional Navier-Stokes equations with mixed derivatives." J. 
Comput. Phys . 41, 1-33 (1981). 

[2] J. S. Allen and S. I. Cheng, "Numerical solutions of the compressible 
Navier-Stokes equations for the laminar near wake in supersonic flow." 
Phys. Fluids 13, 37-52 (1970). 

[3] Yu I. Brailovskaya, "A difference scheme for numerical solution of the 
two-dimensional nonstationary Navier-Stokes equations for a compressible 
gas." Soviet Phys. Dokl . 10, 107-110 (1965). 

[4] L. Crocco, "A suggestion for the numerical solution of the steady 

Navier-Stokes equations." AIAA J . 3, 1824-1832 (1965). 

[5] A. Jameson, W. Schmidt, and E. Turkel, "Numerical solutions of the Euler 
equations by finite volume methods using Runge-Kutta time-stepping 
schemes." AIAA paper 81-1259 (1981). 

[6] R. Peyret and H. Viviand, "Computation of viscous compressible flows 

based on the Navier-Stokes equations." AGARD Report 212 (1975). 

[7] T. H. Pulliam and J. L. Steger, "Recent Improvements in efficiency, 

accuracy, and convergence for implicit approximate factorization 

algorithms," AIAA paper 85-0360 (1985). 



- 21 - 


[8] P. J. Roache, Computational Fluid Dynamics , Herraosa Publishers, 1976. 

[9] R. C. Swanson and E. Turkel, "A multistage time-stepping scheme for the 

Navier-Stokes equations." AIAA paper 85-0035 (1985). 

[10] E. Turkel, "Acceleration to a steady state for the Euler equations." 

Numerical Methods for the Euler Equations of Fluid Dynamics , pp. 281- 
311, SIAM, Philadelphia, 1985. 

[11] E. Turkel, "Algorithms for the Euler and Navier-Stokes equations on 

supercomputers." Progress and Supercomputing in Computational Fluid 
Dynamics (Edited by E. M. Murraan and S. S. Abarbanel), 155-172, 


Birkhauser Publishing Co., Boston, 1985 








Figure lc 



iterations 


Figure Id 







- 29 - 


Figure 


Figure 


Figure 


Figure 


Figure 


Figure 


Figure Captious 

la. Surface pressure distribution for laminar flow over an NACA 0012 
airfoil (M^ = 0.5), Re^ = 5 x 10^, a = 0 degrees). 

lb. Skin-friction (based on free-stream conditions) distribution for 
laminar flow over an NACA 0012 airfoil (M = 0.5, Re = 5 x 10^, 

CO 00 

a = 0 degrees). 

lc. Streamlines for upper surface at the trailing edge (M = 0.5, 

3 

Re = 5 x 10 , a - 0 degrees). 

00 

l d. Convergence histories for laminar airfoil flow calculations. 

A — Runge-Kutta scheme without pseudo-time algorithm (CFL number 
of 1.0). 

B — Runge-Kutta scheme with pseudo-time algorithm (CFL number of 
2.5). 

2a. Surface pressure distribution for turbulent flow over an NACA 0012 
airfoil (M n = 0.5, Re^ = 2.89 x 10^, a = 0 degrees). 

2b. Skin-friction (based on free-stream conditions) distribution for 
turbulent flow over an NACA 0012 airfoil (M = 0.5, 

CO 

Re^ = 2.89 x 10^, a = 0 degrees). 
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Figure 2c. Convergence histories for turbulent airfoil flow calculations. 

A — Runge-Kutta scheme without pseudo-tirae algorithm (CFL number 
of 1.5). 

B — Runge-Kutta scheme with pseudo-time algorithm (CFL number of 


2.7). 
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